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ABSTRACT 

Under ideal MHD conditions the magnetic field strength should be correlated with density in the 
interstellar medium (ISM). However, observations indicate that this correlation is weak. Ambipolar 
diffusion can decrease the flux-to-mass ratio in weakly ionized media; however, it is generally thought 
to be too slow to play a significant role in the ISM except in the densest molecular clouds. Turbulence 
is often invoked in astrophysical problems to increase transport rates above the (very slow) laminar 
values predicted by kinetic theory. We describe a series of numerical experiments addressing the 
problem of turbulent transport of magnetic fields in weakly ionized gases. We show, subject to 
various geometrical and physical restrictions, that turbulence in a weakly ionized medium rapidly 
diffuses the magnetic flux to mass ratio B j p through the buildup of appreciable ion- neutral drifts on 
small scales. These results are applicable to the fieldstrength - density correlation in the ISM, as well 
as the merging of flux systems such as protostar and accretion disk fields or protostellar jets with 
ambient matter, and the vertical transport of galactic magnetic fields. 

Subject headings: diffusion — MHD — turbulence — methods: numerical — ISM:magnetic fields 
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1. INTRODUCTION 

Two dimensionless parameters control the degree to 
which galactic magnetic fields are frozen to the inter- 
stellar gas. One, the magnetic Reynolds number Rm, is 
the ratio of the Ohmic diffusion time to the dynamical 
time, and is typically of order 10^^ - 10^^. The second, 
the ambipolar Reynolds number Rad, is the ratio of the 
ion-neutral drift time to the dynamical time. This num- 
ber is typically many orders of magnitude less than Rm , 
and can approach unity in dense molecular gas. Based 
on these estimates, magnetic fields should be nearly per- 
fectly frozen to the plasma component of the gas, and 
generally quite well frozen to the neutrals, except in the 
densest, nearly neutral regions. 

Thus, the ratio of magnetic fieldstrength to gas den- 
sity B/ p is determined primarily by dynamical rather 
than by microscopic processes. Parameterizing the B — p 
relation by k = din B / din p, one finds k = 1 for com- 
pression perpendicular to B, k = 2/3 for isotropic com- 
pression, K ~ 1/2 for self-gravitating, magnetically sub- 
critical clouds, and k = for compression parallel to B. 

Observations of the B — p relation in molecular gas in- 
deed show that the strong est fields are associated with 
the densest gas ('CmEEifl , iBourke et all l|2001h . 

iSarma et^ l. (2002)). The correlation is consistent with 
K ~ 0.5, although there is so much dispersion, particu- 
larly when upper limits are included, that this relation 
is perhaps only an upper envelope. In atomic gas, the 
B — p relation is consi stent with k ^ over three or- 
ders of magnitude in p ijTroland fc Heileslll98(i() . Obser- 
vational effects alone can introduce substantial scatter 
because all measurements are averages along the line of 
sight and over the telescope beam width, p may not be 
accurately determined, and because only the line-of-sight 
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component of B is measurable. However, even allowing 
for the possibility that the observational scatter is large, 
the B — p relation is strikingly flat. Moreover, there is a 
mean density jump of at least a factor of 50 when going 
from atomic to molecular gas, whereas the correspond- 
ing mean magnetic field strength seems to increase by a 
factor of two or three at most. 

At first sight, the simplest explanation of the observed 
B — p relation is that dense regions in the ISM arise 
primarily from compression parallel to B. Arguments 
against this as the sole explanation are quantitative 
rather than qualitative. In order to coUimate the flow, 
the magnetic energy density should dominate the turbu- 
lent energy density, but the field is at or below equiparti- 
tion. And, if giant molecular clouds are assembled by ID 
compression, they must sweep up material over nearly a 
kpc, too larg e a scale on which to expect coherent flow 
ljMeste nci85V 

In astrophysical environments, the microscopic diffu- 
sivities - whether viscous, resistive, or chemical - are of- 
ten far too small to explain the transport that apparently 
takes place. However, diffusion rates can be enhanced by 
turbulence. Turbulence accelerates transport because it 
creates small scale structure, which diffuses faster than 
the original large scale structure and smoothes the large 
scale gradients. In a turbulent flow with characteristic 
velocity u and correlation time t, the eff ectiv e diffusivity 
Ae is argued to be of order m^t (see eq. jA6p . 

Therefore, in this paper, we investigate whether tur- 
bulence has a similar effect on the transport of magnetic 
fields with respect to the neutral gas. Because the in- 
terstellar magnetic field is subject to ion-neutral drift 
at scales much larger than the resistive scale, we ignore 
resistive effects, except to control numerical diffusion in 
our calculation, and concentrate on the properties of ion- 
neutral drift in a turbulent medium. Although our work 
is motivated by the B ~ p relation observed in the in- 
terstellar medium, it is also relevant to other problems, 
including transport of magnetic fields in weakly ionized 
accretion disks and entrainment of molecular gas by pro- 
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tostellar outflows. 

In we review the theoretical basis for turbulent am- 
bipolar drift and summarize previous results. In ^ we 
introduce our model. Sectionals a description of the nu- 
merical method and its validation. The main results are 
presented in Sj3 Finally, we mention other applications 
in the concluding section, together with a summary. 

2. BACKGROUND 

Consider a medium in which the ionization is suffi- 
ciently low that the neutral and total densities, pn and 
p, are interchangeable, as are the neutral and center of 
mass velocities u„ and u. Under these conditions, the 
continuity equation 



dp „ 



and magnetic induction equation 

dB 

'dt 

can be combined to yield an equation for B/p 



V X (ui X B) 



(1) 



(2) 



|^ + u-v)- = --Vu+ivx(uz5xB), (3) 
dt J p p p 

where u^i = — u„ is the relative drift between the ions 
and the neutrals. The first term on the right hand side 
of equation (PJ represents fieldline stretching, while the 
second term represents ambipolar diffusion. 

In this paper, we restrict ourselves to a geometry in 
which all velocities are in the {x, y) plane, B = zB, and 
all quantities are independent of z. Under these restric- 
tions there is no fieldline stretching, and equation ||3J) 
reduces to 



d_ 

dt 



u- V 



B 
P 



1 



V • [udB). 



(4) 



If the timescales of interest are long compared to the 
ion-neutral collision time , u^i is determined by bal- 
ancing the Lorentz force on the ions against the fric- 
tional force on the neutrals (the so-called strong cou- 
pling approximation). For the geometry assumed here, 
the Lorentz force is due only to the magnetic pressure 
gradient, and 



Substituting equation |(SJ) into equation Q yields 



d \ B 1 

— +u- V - = -V -Xad^B, 

dt J P P 



where 



B' 



Anpiiyi, 



(5) 



(6) 



(7) 



is the ambipolar diffusivity. 

Equation jnj is almost, but not quite, an advection - 
diffusion equation for B/p. The difference is that the 
diffusive flux is proportional to VS, not W{B/p). As an 
immediate consequence, B/p will develop a large disper- 
sion if p its elf is a dvected by the turbulence as a passive 
scalar (see 5.5(1 . 



Equation |(HJ predicts a characteristic ambipolar diffu- 
sion time TAD for a field with characteristic lengthscale 
Lb 



TAD 



^AD 



(8) 



We want to know for which turbulent flow speeds u and 
length scales Lg would such a field diffuse at a rate higher 
than the laminar rate. If the turbulence has a correlation 
time Tc — Lf,/u, the turbulent diffusivity Ag is of order 
u^Tc = LeU (see eq. |A6| ) . The ratio of the ambipolar 
diffusion time tad to the turbulent diffusion time Tt is 
just the ratio of the two diffusivities: tad/ti = Xe/^AD- 
Therefore we are interested in the case Xc/Xad > 1. 
This condition will be satisfied if the magnetic field is 
well frozen to the turbulent eddies. The degree of freez- 
ing is measured by the eddy ambipolar Reynolds num- 
ber RAD{Le, u), the ratio of the ambipolar diffusion time 
across the eddy, TAD{Le) to the eddy correlation time. 
Assuming the eddy size Le is related to u and t^ by 
L„ = uTr we have 



RAD{Le,u) 



Xad 



u 

CAt 



Inserting numerical values. 



RAD{Le,u) = 9.4 X 10-^ u ( ^ 



B J Pn 



(9) 



X,, (10) 



where Lg is expressed in parsecs, u in km s^^, B in Gauss, 
and Un in cm^'^. The p represent molecular weights, and 
Xi is the ionization fraction. The field is frozen to tur- 
bulent eddies which are larger than the size at which 
RAD{Le,u) = 1. For example, in gas with an ionization 
fraction 10^'^, Pi/pn ^ 1, magnetic field B = 5pG, an 
internal velocity dispersion u — 1, and neutral density 
Tin = 50, this critical lengthscale is about lO^'^pc, corre- 
sponding to a column density of about 1.5 x 10^^ cm~^. 

What diffusion rate is actually required to break the 
flux freezing and produce a flat B — p relation? Sup- 
pose a coherent density structure of size L forms on a 
timescale t, with an associated velocity U. The B ~ p 
relation breaks down if the diffusion time is less than the 
formation time. In other words, if the relation 



Xad < LU < Xe 



(11) 



holds, then laminar ambipolar drift is too slow to break 
flux freezing but turbulence is fast enough. 

If L, U, Le, and u are related through the usual scaling 
laws for a turbulent cascade, then U/u > 1 whenever 
L/Lg > 1. Hence turbulent mixing can only destroy 
the B — p relation in large scale structures which are 
controlled by additional physics, such as cooling or self 
gravity, or in regions with strong sources of small scale 
turbulence, such as velocity shear layers. 

Equation © is the basis of three recent studies of 
the effects of t urbulence on the rate of ambipolar drift. 
iZweibell l|2002|) demonstrated accelerated diffusion at ap- 
proximately the eddy rate Ag by finding exact solu- 
tions of equation for random sequences of incom- 
pressible stagnation point flows, valid for a particular 
initial magnetic p rofile and initially constant density. 
iFatuzzo fc Adarna (2002) concentrated on the role of 
fluctuations in B, and hence Xad, in reducing tad- 
They showed that fluctuations in fleldstrength lead to 



Heitsch et al. 



3 



a corresponding; d ispersion in ambipolar diffusion times. 
iKim fc Diamondl (^002) analyzed equation assum- 
ing V • u = and constant p, under the standard as- 
sumptions of quasilinear diffusion theory (QDT, see Ap- 
pendix). They argued that the canonical turbulent diffu- 
sion rate Ae ~ v^t is an upper limit which is approached 
for turbulence well frozen to the eddies. 

Within the framework of QDT, the novel aspect of 
equation (|SJ) is the right hand side, which consists of a 
nonlinear diffusion operator applied to _B, not B j p. The 
latter distinction is unim portant as long as p is uniform. 
IKim fc Diamondl l)2002() argue that diffusion reduces tur- 
bulent transport, because in a highly diffusive system, 
the field is not advected by the flow. While equation ((2Jl 
shows that in the limit \ad = there is no diffusion in 
an absolute sense, we will see that the large scale compo- 
nent can decay by t ransferring powe r to the small scales. 

Taken to gether. iZweibell 12001 . iFatuzzo fc Adams! 
||2P02), and IKim fc Diamondl l|200l support the propo- 
sition that ambipolar drift is accelerated in a turbulent 
medium, provided that the strong coupling approxima- 
tion holds, and in the 2.5D geometry assumed in all three 
papers. The diffusion rate is enhanced by the develop- 
ment of small scale structure in the field, which increases 
the local drift velocity, and by the growth of fluctua- 
tions in fieldstrength, which increases the local diffusiv- 
ity. When the former dominates, the diffusion rate is 
close to the canonical turbulent value u^r, and nearly 
independent of the microscopic diffusivity. 

Although these papers are suggestive, the picture they 
present is still incomplete. Due to the choice of initial 
condition, the difference between diffusion of B and diffu- 
sion oi B/ p\s not addressed by any of the work, although 
the former describes transport of B from an Eulerian vol- 
ume and the latter from a Lagrangian one. The respec- 
tive roles of irreversible flux transport by ion- neutral drift 
and breakup of large scale structure into small scale fluc- 
tuations is not discussed by the QDT calculation, while 
the stagnation point calculation has yet to be embedded 
in a globally valid flow model. In the remainder of this 
paper, we report on a series of numerical experiments of 
the same basic problem which address some, but not all, 
of these limitations. 

3. DESCRIPTION OF THE PROBLEM 

3.1. Assumptions and Equations 

As in fj21 we consider the action of 2D flow (in the (x, y) 
plane) on a perpendicular magnetic field z_B. Because 
we anticipate the formation of structure on small scales, 
we do not assume strong coupling (eq. 0). Instead, we 
compute Ui by solving the ion equation of motion 

Pi ^^Ui + (Ui • V)Ui^ = -^VB^ - PiVin{vLi - u„). 

(12) 

We assume that pi is related to p„ through the 
condition of ionization equilibrium. Although this 
breaks down on short timescales, we have argued else- 
where that it is general ly an excellent approximation 
()Heitsch fc Zweibe]|l2003|) . We also neglect ion pressure 
relative to magnetic pressure. This is a good approxima- 
tion in weakly ionized interstellar gas except near mag- 
netic nulls, which are precluded by the choice of initial 
conditions (see H3 3.5|l . 



We simplify the problem further by treating u„ as 
given. Although this so-called kinematic approach is fre- 
quently made in turbulent mixing problems, it is not 
particularly accurate in the interstellar medium, because 
the magnetic and turbulent pressures are comparable. 
However, in weakly ionized gas there is a lengthscale be- 
low which the neutrals and ions are not well coupled (see 
H3 3.2|l . and one can regard the dynamics as taking place 
below this critical scale. One consequence of treating the 
neutrals kinematically is that p itself behaves as a pas- 
sive scalar. We could reach the same result by assuming 
zero temperature, and neglecting ion pressure. 

With VLi determined from equation (|12|l . we require 
only the magnetic induction equation to close the system. 
In the geometry considered here, the induction equation 
equation (|2Jl is 

= -V • (u,S) + AfjV^B, (13) 

where the second term on the RHS has the same form 
as Ohmic diffusion. We tune this last term so that it 
dominates the numerical diffusion. 

3.2. Linear Theory 

Equations (|12|l and H13|) can be linearized to describe 
small perturbations about a uniform equilibrium state. 
The analytical solutions provide physical insight and 
are useful for numerical tests. Assuming the pertur- 
bations are Fourier modes which depend on (x,y,t) as 
exp [i{kxX + kyy — ut)] yields the dispersion relation 



± 



( — ) 
^ 2 ^ 



(14) 



ky and CAi 



B/^/ Air Pi is the Alfven 



where fc^ = k^ 
speed in the plasma. 

If kcAi ^ I'm/2, the roots of equation H14|l describe 
forward and backward propagating magnetosonic waves 
damped by ion-neutral friction 

Lu X ±kcAi - ■ (15) 

If kcAi < i^in/'^, w is purely imaginary. In the limit 

kcAi ^ Vinjli the roots of equation H14|l are given ap- 
proximately by 

CJ = -^^'^„ (16) 



Ai 



(17) 



Equations H15|l and (jKil) can be compared to the dis- 
persion r elation which accounts fo r ion feedback on the 
neutrals ijKulsrud fc Pearc^ Il969() . According to this 
more complete treatment, w is purely imaginary for 
Si^iTi \/ Pil P < kcAi < Vin/2. The upper limit is exactly 
what is predicted by equation l|14|) . and equation H15|l 
agrees well with the exact dispersion relation. The lower 
limit is not predicted by equation ((111) . At k values well 
below this lower cutoff, the wave propagates at the bulk 
Alfven speed B / ^Airp. In this regime, the neutrals have 
time to be accelerated by the ions within one wave pe- 
riod. This does not happen in our model because we 
omit the drag force on the neutrals. 

The solid line in Figure ^ shows the predicted depen- 
dence of the damping rate (i.e. the imaginary part of 
the wave frequency lo of eq. jl4| ') on Vi„ . The damping 
mode changes at 2cAik as expected. A discussion of the 
overplotted numerical results is left to H44.ll 
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Fig. 1. — Dec ay r ate uji against collision frequency Uin, accord- 
ing to equation I14i for the ID case (solid line) and the 2D case 
(dashed line). The symbols denote numerical results at resolutions 
as labeled in the plot. Note that the 2D case is shifted by exactly a 
factor of 2 with respect to the ID result. The flat parts at low and 
high fin are caused by numerical diffusion, leaving roughly two 
orders of magnitude in in the declining branch of Ui (eq. 1171 ') 
for numerical studies. 



3.3. The Neutral Flow 

We adopt a two-dimensional vers ion of the divergence 
free, "circularly polarized" flow of iGallowav fc Proctod 
l|1992t) . (hereafter GP-flow), which we write in compo- 
nent form as 



—Vf cos(27r(fc/y 4- e sm{2TTkfVft))) 



(18) 



We have chose n this flow be c ause i t has stretching prop- 
erties ( Gallow av fc Proctod l|1992(l show it to be a dy- 
namo flow) and can be written in closed form. Its de- 
gree of chaos can be tuned through the choice of the 
parameter e.^ For e = the flow is steady and possesses 
cellular structure at a single spatial scale. The corners 
of the cells are hyperbolic stagnation points near which 
the fluid undergoes exponential expansion along one axis 
and exponential compression along the other. Fluid cir- 
culates steadily around the center of its cell with an eddy 
turnover time of order 



(19) 



Since each eddy retains its amount of fluid, we expect 
turbulent transport to be minimal. 

The flow pattern for e > can be visualized as eddies at 
a single scale travelling in snake-like patterns across the 
domain. The position of each eddy oscillates by e/kf. 
For e <C 1, this is also roughly the size of the chaotic 
region. Advected quantities are not bound to the ed- 
dies, but travel from one to another, leading to turbulent 
transport. 

No te that e contains a factor of 2tt in our definition 
(eq. US), thus being by th at factor larger than the definition in 
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Fig. 2. — Tracer pattern caused by GP-flow after t = 2Tf. At t = 
0, a tracer quantity is deployed into a circular region with diameter 
L/4. The GP cell number Kf = 5, and e = 0.3. Resolution A'^ = 
1601^. 



Figure 121 shows an example of the GP-flow for Kf, the 
number of cells per domain length L, equal to 5 and 
e = 0.3. Ati = 0, we deploy a tracer quantity into a 
circular region with diameter L/4. The resulting tracer 
distribution is shown at t — 2Tf. We discuss the choice 
of e in mrni 

3.4. Parameter Regimes 

There are four diffusivities: the numerical diffusiv- 
ity Xnum, the magnetic diffusivity An deflned in equa- 
tion H13() . the ambipolar diffusivity Xad deflned in equa- 
tion ITJl, and the eddy diffusivity of the GP flow Ae, which 
is 0{vf/kf) (see eq. \A6\ ). Normalizing by Ag and using 
equation (0), we require 



An 



A. 



A, 



kfCAi CAi 



< 1. 



(20) 



IGallowav fc Proctoil il99^ 



If equation H20() is fulflUed then numerical diffusion is 
smaller than all other diffusion, Ohmic diffusion is weaker 
than ambipolar diffusion, as it is in the interstellar 
medium, and eddy diffusion is faster than ambipolar dif- 
fusion, which allows us to test the hypothesis that turbu- 
lence enhances the rate of magnetic field redistribution. 
The last inequality is equivalent to requiring that the 
field is frozen to the eddies: Rad ^ 1 (see eq. |5]). 

We also seek a separation of lengthscales in the prob- 
lem. The largest scale possible in a periodic domain is 
L, while the eddy scale is L/ Kf. The distinction between 
large and small scales can be maintained only if 

Kf » 1. (21) 

Finally, as we saw in equation (|14|l . waves propagate 
only for sufficiently large wavenumbers. We are actually 
interested in the nonpropagation regime, because only in 



Heitsch et al. 



5 




Table 1. Ambipolar Reynolds number Rad 



Fig. 3. — Valid regimes (schematically) for the ISM (horizontal 
lines) and the models (diagonal lines) under the strong coupling 
condition (eq. |5]). Ve rtica l lines denote the nonpropagation limit 
according to equation I22i for the ISM and the models. 



that case is the neutral flow imprinted on the ions. Thus, 
we require 



kfCAi 1 
^ 9- 



(22) 



When equation (|22|l is combined with the last condition 
in equation (|20ll . the result is an upper limit on CAi/vf. 
In the interstellar medium, where the magnetic field is 
close to equipartition with the turbulence, this quantity 
can be of order (p/p.;)^/^, and hence quite large. For 
computational reasons, we are unable to make kfVf/uin 
as small as it should be, and therefore CAi/vf will be 
smaller than it ought to be. However, there is no physical 
reason why turbulence cannot obey the proper ordering. 

The restrictions on Vf/cAi and kf CAi/vin resulting 
from equations (|20|l and (|22|l are sketched in Figure O 
For kf CAi/vin > 1, we are in the weakly damped branch, 
and anything below the diagonal thick line corresponds 
to Rad < 1- 

Taking these constraints together, the range of permis- 
sible i>in is limited to at most 2 orders of magnitude, the 
range of B by slightly more than one order of magnitude, 
and Kf between 5 and 10. 

Tables Hand 12 summarize the parameters for all mod- 
els. For each combination of B and Vin, we ran models 
for three GP forcing scales, namely Kf G {5, 7, 10} at a 
resolution of = 801^. These models were completed by 
a set of pure tracer models (i.e. models with = 0, and 
B subjected to turbulent diffusion only). To assess reso- 
lution effects on the results, we repeated the calculations 
at iV = 1601^. 

3.5. Initial Conditions 

We prepare our system with the initial conditions as 
given in equation H23|l for the 2D case. 



IB Ui„ - 


0.7 


2.3 


7.1 


23.0 


0.1 


3.4 X 10^ 


1.1 X 10'' 


3.4 X 10"* 


1.1 X 10^ 


0.5 


1.3 X 10^ 


4.5 X 10^ 


1.3 X 10*^ 


4.0 X iU 


1.0 


3.4 X 10^ 


1.1 X 10^ 


3.4 X 10^ 


1.1 X 10^ 


2.0 


8.6 X lO'^ 


2.8 X 10^ 


8.6 X 10^ 


2.8 X 10^ 




Table 2. Ratio A 


AD /An 




IB Uin ^ 


0.7 


2.3 


7.1 


23.0 


0.1 


4.8 X 10^1 


1.4 X 10-1 


4.8 X 10-2 


1.4 X 10-2 


0.5 


1.2 X 10^ 


3.6 X 10° 


1.2 X 10'^ 


4.0 X 10-1 


1.0 


4.8 X 10^ 


1.4 X 10^ 


4.8 X 10° 


1.4 X 10° 


2.0 


1.9 X 10^ 


5.8 X 10^ 


1.9 X 10^ 


5.8 X 10° 



B{x, y, 0) = Bo + -Bi [1 -I- cos{kpx) cos(fcpy)] 

p{x, y, 0) = po - ^pi [1 + cos(fcpx) cos{kpy)] . (23) 

The ID case is initialized similarly, with y = 0. 
We add the offset in B to prevent field reversals, as 
these are known to lea d to singularities in the pres- 
ence of ambipolar drift ijBrandenburg fc Zweibell 119951: 
'Mac Low fc Smithlll997D . We choose the perturbation 
amplitudes to be Bi = O.liJo, Pi = O.lpo, and the per- 
turbation wave number kp to be 27r/L, i.e. the pertur- 
bations reside on the largest possible mode. According 
to the formalism of QDT discussed in the Appendix, the 
kp components of B and p represent the large scale fields 
which are predicted to diffuse under the action of the 
small scale fields. 

Ionization equilibrium f ^3 3.1|) implies that pi is slaved 
to p. In the physical environments of interest here, 



Pi (X p 



1/2 



Thus, the 10% perturbation of p imposed 
here (eq. [2S1) would cause a 5% perturbation of pi. For 
the sake of simplicity in the numerical scheme, we kept 
Pi constant in time and space. We believe the 5% vari- 
ation in ion density to be too small to affect our results 
qualitatively. 

To summarize, we initialize B and p using equa- 
tions H23I) and evolve the system in space and time using 
equations H12|) and p3(l . where u„ in equation 112|l is 
given by equation (|18|l . 

4. NUMERICS 

4.1. The scheme 

We solve equations (|12() and (|13|l using a modified ver- 
sion of the 1st order g as-ki netic flux-sp l itting method de- 
scribed by|x3 1(19991) and iTang fc XTil MUSCL^ 
limiters allow a 2nd order reconstruction of the flow 
variables at the cell walls. As we are evolving only 
zB{x,y,t), V • B = is satisfied. The (otherwise pas- 
sive) z velocity component serves as a tracer in order to 
measure the transport properties of the flow. The AD 

^ Monotone Upwind Schemes for Scalar Conservation Laws 
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drag force terms are added as external source terms to 
the cell-centered momenta a nd can be treated implic- 
itly if necessary (similarly to lMac Low fc SmithI l)1997D : 
|T6th (1995)). The CFL-condition includes the diffusion 
timestep. As we are interested in representing the peak 
value of the initial magnetic field exactly, we employ a 
grid with an odd number of cells (the symmetry axis of 
the initial magnetic field goes through a cell center). 

Figure ^ serves as a validation of the method for our 
problem. The decay of the magnetic field is measured in 
terms of the peak magnetic field amplitude against time. 
The initial conditions for ID and 2D are given in equa- 
tions (|23|l . Apart from the fact that the numerical data 
points agree with the predictions in most cases at the 2%- 
level, we note the flattening at low and high Vin. Both 
are caused by numerical diffusion. For small z^i„, this dif- 
fusion arises from the limited grid resolution. At Vnum, 
the point where the numerical results start to leave the 
predicted curve, the physical diffusion becomes smaller 
than the intrinsic numerical diffusion, which from then 
on governs the diffusive behavior. As soon as the numer- 
ical curve flattens, the results are meaningless for any 
physical problem. The timestep in this weakly damped 
branch is controlled b y the diffusion timestep as in e.g. 
iMac Low et"all lll995j) . 

For large the diffusion is caused by the finite 
timesteps. If the collision time becomes smaller than 
the timestep At, the number of collisions per At is lim- 
ited. In other words, At <C must be guaranteed in 
order to prevent losing collisions and thus introducing nu- 
merical diffusion. The timesteps would be impractically 
small, and pursuing the computations would require an 
implicit treatment of the whole ion momentum equation. 
Although this is theoretically possible, we chose to im- 
plement a 2nd order Runge-Kutta time stepping method 
which guarantees results unaffected by numerical diffu- 
sion up to Vin ~ 70 (see Fig. P). 

Note that the predictions for the 2D decay rate (dashed 
line) coincides exactly with the numerical datapoints as 
well. It is this agreement which makes us confident in us- 
ing the method for the investigation reported here. Fur- 
ther validation is presented in fjs] 

4.2. Control of numerical diffusion 

As we mentioned in connection with equation , we 
control the diffusive behavior of the scheme at grid scale 
by introducing an artificial resistivity An which dissipates 
the field above the grid scale. As is usually the case 
in astrophysical computations, Ao is much larger than 
the physical resistivity. However, it is necessary because 
the GP-flow mixes scales very efficiently, leading to dif- 
fusion at the smallest possible scale after approximately 
an eddy-turnover time r/ . This artificial resistivity guar- 
antees well-behaved diffusion properties above the grid 
scale. It acts both on the magnetic field and the tracer 
field, which we identify as p. Thus, B and p are guaran- 
teed to be smoothed on the smallest scales in the same 
manner. 

In the flux-splitting scheme as described above, the 
amount of diffusion (last term of eq. is given by the 
slope of -B, which in turn is computed during the recon- 
struction of the magnetic field at the cell walls. Thus we 
can implement the resistivity efficiently in our scheme. 
Restricting ourselves for the moment to a purely resistive 
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Fig. 4. — Dec ay rate LOi against Ohmic diflfusivity Aq according 
to equation I24i (solid line) for the ID case and (dashed hne) for 
the 2D case. The symbols denote numerical results at resolutions 
as labeled in the plot. Note that the 2D case is exactly a factor of 
2 shifted against the ID result. The flat parts at low An are caused 
by numerical diffusion. 



medium - excluding ambipolar diffusion - the dispersion 
relation in the physically interesting regime (small An) is 
given in ID by 



uo = ± 



'Ai 



fc2 - ( 



Aofc2 , 



(24) 



where we are interested in the case \nk <C cji^fc • The 
corresponding tests are shown in Figure ^ Numerical 
diffusion (flat part for small An) occurs at the same level 
as for the implementation of ambipolar diffusion (see 
Fig. Comparing the values of oji in Figure ^ with 
those of Figure n shows that our parameter space is lim- 
ited by resistive diffusion. 

5. RESULTS 

We computed the evolution of B and of a tracer field for 
a set of realizations of the basic model described in fj^l see 
Tables 1 and 2. We found that some, but not all, aspects 
of the transport of these quantities can be cast in terms 
of turbulent diffusion. After explaining how we measure 
the total diffusivity D in our models (see Appendix), we 
discuss the transport properties of the GP-ffow in MS 5.11 
In H5 5.2l and ii5 5. 31 we describe the evolution of the tracer 
p and the transport of B, respectively. Turbulent trans- 
port affects p and B in different ways, requiring a closer 
view on the ion ffow (|5 5.41) . Finally in il5 5.5l we discuss 
the evolution of the flux to mass ratio B/ p. 

5.1. Transport Properties of the GP-flow 

The transport properties of a flow are closely re- 
lated to its Lyapunov exponents (see for example 
jQiaziiJ 1122^). We used a meth od impleme nted by 



Brumell. Cattaneo fc Tobia^ l)200lD . based on iSowardI 
l)1994|) to calculate the Lyapunov exponents for the GP- 
flow. It follows the time evolution of the separation d of 
two points, initially an inflnitesimal distance do apart. If 
one or both points are located in a chaotic region, the 
separation grows exponentially. The grid resolution (in 
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Fig. 5. — Mean and maximum Lyapunov exponents against e. 
Error bars show the standard deviation across one frame. The 
maximum exponent Amax peaks at small e because in this case, 
particle orbits are long and the chaotic regions are located at the 
cell boundaries. The flow's transport properties are determined 
primarily by the mean exponent (A). 



this case 128^ and 256^) is given by the number of paired 
starting points for foUowing the particle trajectories, in 
our case distributed evenly over the flow domain. 

Figure El gives the dependence of the domain averaged 
Lyapunov exponent (A) and maximum Lyapunov expo- 
nent Amax on e. The mean Lyapunov exponent (A) peaks 
near e « 0.25, corresponding to a maximum phase shift 
of 7r/2 (see eq. llSj). The maximum exponent shows con- 
siderable scatter (mirrored in the large dispersions about 
the mean), but it too drops nearly monontonically for 
e > 0.25. This came somewhat as a surprise to us, as a 
rising e increases the time-dependence of the GP-flow. 

Figure El explains this behavior. For e > 0.25 (maxi- 
mum phase shift tt/2), the frequency in velocity pertur- 
bations is doubled, disrupting the coherent flow struc- 
tures and reducing the distance a test particle can be 
transported. Based on this qualitative argument, we ex- 
pect a maximum in the transport rate for e ~ 0.25. 

Figure [3 demonstrates the connection between the 
complexity parameter e in the GP-flow (eq. ^Hl) and the 
diffusion constant D (eq. |A9p . The diffusion constant 
peaks approximately at the e for which the mean Lya- 
punov exponent (A) is maximal (see Fig. jSJ . The error 
bars reflect temporal fluctuations in D. Because we wish 
to maximize the transport, we chose e — 0.3 in all the 
models discussed subsequently. 

We assessed the effect of the GP-flow's temporal peri- 
odicity on the transport rate by running some numerical 
experiments in which we replaced the factors of kfVft in 



Fig. 6. — Velocities Ux and Uy according to equation I18i at 
X = y = 0. For e > 0.25 (maximum phase shift tt/2), oscillations 
in Uy change to the next higher octave, the same happens for n^, 
for e > 0.5. By this, the coherent flow structure is disrupted. 
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FlG^ 7. — Diffusion constant D against complexity parameter e 
(eq. 1181 ^. The er ror b ars give the standard deviations about the 
mean in equation IA9i . 



equations l|18|) by random phases chosen from a uniform 
distribution in (0, 27r) at fixed intervals in time. The re- 
sulting transport was less than for the GP flow with the 
same value of e, so we pursued this model no further. 

5.2. Transport of a Tracer 

The transport of a passive scalar by the GP flow es- 
tablishes a baseline against which the transport of the 
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Fig. 8. — (k = l)-modc of Q - where Q is either the density p or the magnetic field strength B - against time for resolution A'^ = 801^, 
K/ 6 {5, 10} and e = 0.3. (a) Decay of density field. The solid {nf =5) and dashed [kj = 10) lines stand for the turbulent decay rate 
expected from QDT (eq. IA6I V (b) Decay of magnetic field with = 0.7. (c) Decay of magnetic field with Ui„ = 23.0. The slope gives the 
decay rate uji. Panels (d) to (f) show the difference (Qgoi — Qi60i)/Qi60l in percent for runs of resolutions = 801'^ and N = 1601^ and 
parameters corresponding to panels (a) through (c). For Ui„ = 0.7, B develops strong peaks, which tend to be underresolved at N = 801'^. 



magnetic field can be compared, ft is also the relevant 
transport equation for the gas density p. We compute the 
transport by integrating equation (|f 3|l with = u^p. 
The initial conditions are given by equations H23|) . Nu- 
merical diffusion is controlled as discussed in ij44.2l 

Panel (a) of Figure |H1 shows the time evolution of the 
average of the F ourier coefficients of pk=i (as defined fol- 
lowing eq. \A8\ ) for two GP cell sizes (k/ = 5, 10). The 
curves clearly demonstrate exponential decay. The differ- 
ence in decay rates between the Hf — 5 and kj — 10 cases 
agrees reasonably well with the expected u^t scaling of 
the turbulent diffusivity, since doubling Kf decreases r 
by a factor of 2 (eq. ^HI)- We note that the deca y ra tes 
do not reach the values predicted by QDT (eq. |A6p ''. 
The bumps in the curves are caused by time variations 
in the GP flow. Not only does the temporal frequency 
of the bumps correspond to the GP frequency, but the 
difference between two spatial resolutions (Fig. |SId)) is 
very small over the duration of the calculation. 

Figure El gives D as defined in equation (|A9|) for the 
same GP tracer models as in Figure|Hl The settling effect 
of larger GP cell numbers is clearly visible. Note again 
that the different resolutions give identical results. 

The top row of panels in Figure lTUI show Fourier spectra 
of p for £ {5, 7, 10}. The vertical dashed lines indicate 
multiples of Kf. All spectra display at least two branches 
with significantly different power. The values in the max- 
imum branch group at riKf ± 1, with n S {1, 2, ...}. This 
is a direct consequence of the nature of the advection 
operator. The velocity u with its single scale Kf initially 
beats together with the density perturbation, which is 

This does not come as a surprise since QDT provides an es- 
timate whose general assumptions are quite different from the as- 
sumptions made here: e.g, the GP-flow has a single scale and it 
is spatially and temporally periodic. The oscillation of the eddies 
in position prevents the maximum stretching of fluid elements in 
GP-cell corners. 
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Fig. 9. — Diffusion constant (eq. IA9I ) for the turbulent decay of 
the passive scalar, identified with the density p. The open squares 
denote results for turbulent ambipolar diffusion at the highest col- 
lision frequency (ui„ = 23.0), demonstrating that for large collision 
frequencies the ions follow the GP-flow closely. 



at A: = 1 (see eq. This creates density structure 

at Kf ± 1. Advection of this secondary structure gen- 
erates power at 2Ky ± 1, and so on. The spectra should 
have power at these wavenumbers only; the power seen in 
Figure El at other k (all branches except the maximum 
branch) is entirely due to numerical noise*. The ampli- 
tude of the noise is 2-3 orders of magnitude less than the 
power at wavenumbers UKf ± 1 (the maximum branch) 
up to the Nyquist fc, which we regard as acceptable. 

Although eqs. IT2t and CH} are solved in double precision, the 
Fourier transform is computed in single precision. 
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Fig. 10. — Fourier spectra of p (upper row) and B (lower rows for minimum and maximum collision frequency Vi^) for k/ G {5, 7, 10}. 
Vertical dashed lines denote multiples of kj (see text). Resolution N = 801^. Each of the nine plots contains one spectrum. The Nyquist 
fcjv = 400. 



5.3. Transport of the magnetic field 

Panels (b) and (c) of Figure |S1 show the decay of the 
(±1, ±1) Fourier coefhcients of B for two resolutions and 
two collision frequencies. The decay rates, when aver- 
aged over the GP flow period, are well fit by exponen- 
tials, and the scaling of the rates with k/ is consistent 
with turbulent diffusion. The decay rate at the higher 
collision frequency is close to the GP rate shown in the 
leftmost panel. At the lower collision frequency, the field 
decays substantially faster because of the relatively large 
value of Xad in that case: the turbulent and laminar 
diffusivities are additive. The difference between these 
cases reflects the greater degree to which the ion flow is 
slaved to the neutral flow in the high collision frequency 
case. 

Panels (e) and (f) of Figure |H1 give the differences 
between these calculations for two resolutions. At the 
higher collision frequency, the differences increase grad- 
ually with time and remain small over the duration of 
the calculation (although they are larger than the corre- 
sponding differences for the tracer p shown in the lower 
left panel). At the low collision frequency, the difference 
grows quickly and saturates at about 10%. This shows 
the increased role of small scale structure in at low 
collision frequencies. Smaller lead to larger compres- 
sion and thus stronger peaks in B, which in turn tend to 
be underresolved at = 801^. 

The diffusivities D computed according to equa- 



tion ljA9|l for the two resolutions are nearly indistinguish- 
able (Fig. The good agreement between the theo- 
retically predicted quiescent diffusion rates Xad and the 
numerical result validates our method of measuring D as 
well as provides an additional test of the underlying nu- 
merics. The approximate invariance of D over time, and 
insensitivity of D to show that the transport of B 
is well described as turbulent diffusion (the small offset 
between the turbulent diffusivity for = 0.7 and the 
other cases is again caused by the relatively large value 
of Xad in that case). 

Because the ion flow is modified by Lorentz forces, 
we investigated the dependence of D on magnetic field- 
strength. FigurclT^plots the time-averaged D - denoted 
(D) - for four collision strengths and four magnetic field 
strengths. The laminar diffusivities (lines) follow the 
expected dependence of diffusivity on B. In contrast, 
the turbulent diffusivities (calculated by measuring D in 
a turbulent model and subtracting Xad) depend only 
marginally on B, despite the fact that the Alfven Mach 
number varies from less than a tenth to slightly greater 
than unity. Note that at the lowest value of B, the quies- 
cent diffusion rates are nearly independent of at this 
value of B, the quiescent diffusivities are completely con- 
trolled by An- On the other hand, at the highest value of 
B and lowest value of Vin, the laminar diffusivity domi- 
nates the turbulent one. This case is less interesting to 
us for studies of the ISM. 

The Fourier spectra of B are shown in the second and 
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third rows of Figure E| for two collision frequencies and 
K/ G {5,7,10}. In striking contrast to the spectra of 
p (top row of this figure), the spectra peak at multi- 
ples of K/. This effect is particularly strong at the low 
collision frequency, and provides graphical evidence that 
the transport of B to small scales is quite different from 
that of p. The disparity must be caused by differences 
between and the GP flow, i.e. the ion-neutral drift. 
To see this, we rewrite the induction equation H13|) . with 
resistivity included, in the form 



AoV^ B = -BV • u.,: 



(25) 



The second term on the LHS of equation (|25l) represents 
advection. It is predominantly, although not entirely, 
advection by the GP flow. As such, it cascades power 
to TiKf zL 1, just as occurs for the tracer. The RHS of 
equation H25|l represents compression. Here the GP flow 
plays no role. If the inhomogeneous part of B were less 
than the mean part, the RHS could be approximated 
by BqW ■ Ui and the Fourier spectrum of V • U; would 
map directly to that of B. Since B is roughly twice Bq, 
this is only part of the story; the compression term gives 
significant nonlinear coupling. In any case, it is clear 
that we must investigate the ion flow. 

5.4. The Ion Flow 

The ion flow is driven by friction with the neutrals, 
by magnetic pressure, and by its own Reynolds stress. 
Friction drives toward ugp- The Reynolds stress is 
nonlinear, and drives compressive flow at multiples of 
Kf. The magnetic field begins with power in the (1,1) 
and (0, 0) components. This power is cascaded to smaller 
scales by advection, but is also coupled to the compres- 
sive part of the ion flow. The Lorentz force itself is non- 
linear, creating additional harmonics in the ion flow. 

We measure the compressibility of the flow through the 
parameter (i?), the ratio of the rms divergence of to 
its rms curl. Figure IT^ shows (R) as a function of for 
four magnetic fleldstrengths. At the lowest collision 
frequency the divergence is a few tenths the curl, but 
quickly decreases with increasing i^in- The value of B 
influences (R) less than the value of i^m does, and (R) 
is not monotonic in B. At low B, Lorentz forces are 
insigniflcant and compression arises from the Reynolds 
stress. As B increases, Lorentz forces begin to play a role 
and lead to some additional compression. As B increases 
further, the field resists compression and actually reduces 
Vui. 

The spectra of the u^^; and V • are shown in Fig- 
ures El and El • There is substantial power in both 
quantities at k/ and its multiples. This can be attributed 
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primarily to frictional driving by uqp with generation of 
power at the multiples of k / by Ui ■ Vui . The peak in the 

spectrum of Uix at Kf reflects the close correspondence 
between the ion flow and the GP flow. The spectrum of 
V • Ui, however, is peaked not at Kf, but at its second or 
fourth multiple. The lower the collisionality, the higher 



the peak. The flow generated by the nonlinear Reynolds 
stress as the GP flow beats together with itself is com- 
pressive, and these flows themselves generate higher har- 
monics yet. If there were no Lorentz force, the spectrum 
of Uj would consist purely of harmonics of the GP flow. 
The magnetic fleld, however, is cascaded by nonlinear 
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Fig. 16. — Magnetic field strength B against density p after 
one eddy turnover time Tf for kf = 5 and Ui„ = 0.7. The line 
denotes the initial condition, and simultaneously the result when 
B is subjected to turbulent diffusion only, but not AD, after one 



advection to the sideband wavenumbers riKf ± 1. The 
resulting power in the magnetic pressure gradient fills in 
the spectrum of at other wavenumbers. 

In the interstellar medium, the respective magnitudes 
of the Reynolds stress and Lorentz force terms may well 
be reversed. As we commented below equation (j22|l . nu- 
merical considerations compel us to make Vf /cAi unreal- 
istically large. However, both these terms are nonlinear, 
and both cascade the ion flow to small scales. 

5.5. Diffusion oj B / p 

Finally, we turn to the B — p relation itself. A direct 
demonstration of how quickly B / p changes for the GP- 
flow in combination with AD is given in Figure lTEl After 
one eddy turnover time, the initially strict correlation 
between B and p (solid line) is completely destroyed. 
The solid line also represents the correlation B{p) for B 
subject to the GP-flow alone, not including AD. In that 
case, B and p are both tracer fields, and diffuse in the 
same manner. 

The separation of B and p on small scales is reflected 
in the decay of the largescale, or (±1, ±1) component of 
B/p (Fig. I17|) . Since - neglecting numerical diffusion - 
p ranges only from 0.9 to 1.1, we have to a fairly good 
approximation 

\P J k=l Po Po Po 

Since both Bk=i and pk=i decay (see Fig. ITTjl . it is 
clear that {B/p)k=i must decay as well, at least in the 
limit represented by equation l|26|) . The decay rate esti- 
mated by equation H26|l approximates the actual decay 
rate quite well. 

6. SUMMARY 

The magnetic fieldstrength - density relation in the in- 
terstellar medium is observed to be quite weak, particu- 
larly in the atomic component and through the transition 
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Fig. 17. — Decay of the (k = l)-component of B and B/p against 
time for Uin G {0.7, 2.3}. 



from diffuse gas to GMCs. The flatness of the trend is 
consistent with field-aligned flow or with an enhanced 
diffusion rate. This paper is concerned with testing the 
hypothesis that the diffusion rate is due to ambipolar 
drift accelerated by turbulence. We report on a series of 
numerical experiments intended to complement, and in 
some ways extend, analytical work on this problem by 
Fatuzzo fc Ada ms (2002), Kim & Diamond ( 2002), and 
IZweibell ([20021) . Within the framework of interstellar 
medium physics, the setup is best imagined as small scale 
turbulence acting on a relatively well ordered field. 

We assumed a 2.5D geometry in which the motions 
are perpendicular to the magnetic field zB and indepen- 
dent of z. In this situation, the fieldlines are shuffled but 
not bent, and there is no magnetic tension force. This 
type of turbule nce is expected in a strong magnetic field 
JStrausf^'igTe*) or as the outcome of an MHD cascade 
(Goldrcich & Sridhar 1995), although in completely sup- 
pressing variation with z we have taken this description 
to an extreme degree. 

We imposed a f low o f the type considered by 
iGallowav fc Proctoill)1992|) (see eq. ^]) on the neutrals. 
The GP flow is spatially and temporally periodic, diver- 
gence free, and chaotic. This flow lacks some obvious fea- 
tures of interstellar turbulence, notably compressibility 
and a full spectrum of spatial scales. We chose this flow 
because it can be written in closed form, its level of chaos 
can be tuned by specification of a single parameter, and it 
is known to produce eddy diffusion. This made it a good 
candidate for a first series of numerical experiments. The 
basic assumption that the neutral flow can be prescribed 
as independent of Lorentz forces can be justified in either 
one of two ways. If the magnetic field were well below 
equipartition with the turbulent flow, this assumption 
would be valid independent of scale. This is, however, 
generally not the situation in the interstellar medium. 
Alternatively, we could restrict the computation to eddy 
sizes less than the cutoff wavenumber for strongly cou- 
pled MHD wave propagation; k > 2{pi/ pnY^^iym/cAi 
We emphasize, however, that this restriction is only a 
self consistency requirement for our calculation, not for 
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turbulent ambipolar diffusion itself. 

We assumed a weakly ionized gas in which the bulk 
density p is advected with the flow equal to the GP flow 
and the ion density pi is kept constant (see i j3 3.5|l . We 
solved the ion momentum equation (eq. |12p including 
the full nonlinear advection operator, the magnetic pres- 
sure gradient, and friction with the neutrals. This al- 
lowed us to go beyond the strong coupling approximation 
(eq. [S]) made in the analytical work. However, we chose 
the GP wavenumber k f low enough that there would still 
be substantial friction in an eddy turnover time (eq. \22\). 
The magnetic field is advected and compressed by the ion 
flow, and subject to a small amount of resistive diffusion, 
according to equation itT^ . 

Our primary consideration, in choosing the computa- 
tional scheme, was to minimize the role of numerical 
diffusion. Tests of the code based on linear theory, on 
measured rates of laminar diffusion, and with regards to 
resolution are shown in Figures |H1 and By scrupu- 
lous choice of parameter regime we are able to keep the 
anomalous growth of small scale structure to a few per- 
cent or less while still considering collision frequencies 
varying by a factor of slightly more than 30, magnetic 
fieldstrengths over a factor of 20, and up to a factor of 
10 between the large scale quantities to be diffused and 
the GP eddies. 

The most serious physical shortcomings of the model, 
as a realization of interstellar turbulence, are its re- 
stricted geometry, kinematic prescription for the neutral 
flow, and relatively large ratio of turbulent speed to ion 
Alfven speed. Although all three of these features were 
to some extent forced upon us by numerical consider- 
ations, we see some advantages to the first two, which 
have enabled us to study turbulent ambipolar diffusion 
in a relatively uncomplicated setting without a bevy of 
competing physical effects. The main consequence of the 
third feature is that it makes the Reynolds stress overly 
prominent in comparison to the Lorentz stress. However, 
even when V • is relatively small, the differences be- 
tween Ui and ugp are significant enough on small scales 
to decorrelate B and p (Fig. EJ- 

These are the main results of our calculations: 

1. The GP flow causes decay of the large scale compo- 
nent of the density (or any other tracer field) . The 
decay is brought about by mixing to small scales. 
Despite the relatively small number of eddies in the 
system, the decay rate is well described by a tur- 
bulent diffusivity of order Vf/kf, as expected from 
mixing length theory. The decay of the tracer is 
shown in Figure |S1 

2. The ion flow driven by the GP flow causes decay of 
the largescale component oi B. At the largest value 
of Vin in our experiment, the decay rate is close to 
the decay rate of the tracer. At the smallest value, 
it is faster. We attribute this, at least in part, 
to the substantial role of laminar diffusion in this 
case. The eddy diffusivity is nearly independent of 
the microscopic diffusivities (Ohmic and ambipo- 



lar) , and of the numerical resolution. The decay of 
the fleld is shown in Figure |H1 

3. Although the large scale structure in p and B decay 
at similar rates, the large scale structure in their ra- 
tio, B/p, decays at the eddy rate as well, as shown 
in Figure El This is brought about by the differ- 
ences between ugp and on small scales. These 
small scale relative drifts rapidly destroy any cor- 
relation between B and p, as shown in Figure [TCI 
Thus, neither point to point measurements nor line 
of sight averages would yield a B — p relation. We 
regard the separate transport of B and p as Eule- 
rian transport, and the decay of B/p as Lagrangian 
transport. 

What are the implications for the interstellar medium? 
The enhanced diffusion rate demonstrated here must be 
balanced against large scale compressive flows, which act 
to restore the B — p relation. There is evidence that in 
environments such as H I shells, which are produced by 
strong dynamical compression and which show relatively 
strong magnetic fields, turbulence is secondary to flow. 
This may also be the case in dense molecular gas, in 
which frictional damping of the turbulence is too strong 
to permit the small scale ion-neutral drifts necessary for 
diffusion. In future work, we intend to include full neu- 
tral dynamics, which would allow such flows, driven, for 
example, by cooling, or by self gravity. This would also 
allow an improved realization of the turbulent spectrum. 
Inclusion of a third dimension would permit us to exam- 
ine the role of field aligned flow, and to consider stretch- 
ing as well as compression of the field. Both can affect 
the B — p relation. 

Finally, we mention some other applications of turbu- 
lent diffusion of the magnetic field with respect to the 
gas. It may play an important role in the escape of the 
large scale horizontal field from the Galactic disk. It may 
also permit the mixing of stellar fields with in situ fields 
in weakly ionized accretion disks, and jet fields with am- 
bient fields in outflows from young stars. 
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APPENDIX 

MEASUREMENT OF TURBULENT DIFFUSION 

Most discussions of turbulent diffusion are based on quasilinear the ory. We provid e a brief review of quasilinear 
diffusion theory here; for details and a more rigorous derivation see e.g. Moffatt (TmB)- 
Consider a quantity q which evolves according to the advection - diffusion equation 

^ + u • g = V • AVg (Al) 

where A is the microscopic diffusivity. Assume q can be decomposed into a mean part (q) , the ensembl e ave rage of q, 
and a fluctuating part 6q with zero mean, while u has zero mean and is isotropic. Averaging equation IjAip gives 

^ = -(u.V<5q)+V.AV(<z). (A2) 

We solve for 6q by subtracting equation ljA2|l from equation (|A1(I and discarding the terms quadratic in the fluctuations. 
The result is 



^ = -u- V((7> + V • AVfJq. (A3) 
ot 



If the correlation time r of the turbulence is short compared to the diffusion time, the solution of equation (|A3|) is 
approximately 



Sq 



J u-V{q)dt'. (A4) 



Substituting equation HA4|I into equation HA2|I and averaging over an ensemble gives 

d{q) 



dt 

where 

,2 



V- (Ae + A)V(<7), (A5) 
Ae - uV (A6) 



is the turbul ent diffusivity. 

Equation (|A2|) shows that turbulent transport is the result of the turbulent velocity field u beating together with 
fluctuations in the field to be transported, in this case q. With the assumption that the scales of the mean and 
fluctuating fields are well separated, Sq is proportional to the local gradient of (q), and the turbulent fiux has the 
same form as a diffusive flux. An immediate implication of equation (|A2|I is that if {q) is initially a Fourier mode; 
(q) (x, 0) ~ (q) exp (ik • x) , then (q) decays exponentially, with 

(9)(x,t) = (q)exp(-Ae/c^i + ik-x). (A7) 

In order to in vestigate whether large scale quantities undergo turbulent diffusion in our simulations, we integrate 
equation (jA5|l over a domain of area A with boundary contour C. Using Gauss's theorem we derive 



d_ 

dt 



■ / {qu=i)dxdy = (Ae + A) / n • V{qk=i)dl, (A8) 

J A Jc 

where we assu me Ag and A are spatially constant. The total diffusion coefficient D is the ratio of the left hand side of 
equation (jA8|l to the right hand side. 

We isolate the large scale fields in our model by setting all modes in the Fourier transform to zero, except the ones 
corresponding to kx^y = ±1. We label the resulting quantity qk=i- We identify the average of qk=i over one GP 
flow period t^, {qk=i), with (q). Since the domain average of (q) is zero, we choose A to cover one quarter of the 
computational domain, cente red on the maximum of the initial magcntic field perturbation. Computing the left and 
right hand sides of equation (|A8p leads to an expression for D 

D=j^ l^{qk=i)dxdy/ J^n-V{qk=i)dl. (A9) 



REFERENCES 



Bourke, T. L., Myers, P. C, Robinson, G., & Hyland, A. R. 2001, 

ApJ, 554, 916 
Brandenburg, A. & Zweibel, E. G. 1995, ApJ, 448, 734 
Brummell, N. H., Cattaneo, F., Tobias, S. M. 2001, Fl. Dyn. Res., 

28, 237 

Crutcher, R. M. 1999, ApJ, 520, 706 

Drazin, P. G. 1992, Nonlinear Systems, pp 140 - 143 (Cambridge: 
Cambridge University Press) 



Fatuzzo, M. & Adams, F. C. 2002, ApJ, 570, 210 
Galloway D. J. & Proctor, M. R. E. 1992, Nature, 356, 691 
Goldreich, P.M. & Sridhar, S. 1995, ApJ, 438, 763 
Heitsch, F. & Zweibel, E. G. 2003, ApJ, 583, 229 
Kim, E.-J. & Diamond, P. H. 2002, ApJ, 578, 113 
Kulsrud, R. & Pearce, W. P. 1969, ApJ, 156, 445 
Mac Low, M.-M, Norman, M. L., Konigl, A. & Wardle, M. 1995, 
ApJ, 442, 726 



Heitsch et al. 



Mac Low, M.-M & Smith, M. D. 1997, ApJ, 491, 596 

Mestel, L. 1985, in Protostaxs and Planets II, ed. D. C. Black & 

M. S. Matthews (Tucson: Univ. Arizona Press), 320 
Moffatt, H.K. 1978, in Magnetic field generation in electrically 

conducting fluids, Cambridge Univ. Press, Cambridge 
Sarma, A.P., Troland, T.H., Crutcher, R.M. & Roberts, D.A. 2002, 

ApJ, 580, 928 

Sowaxd, A. M. 1994, in Lectures on Solar and Planetary Dynamos, 
Publications of the Newton Institute, Vol. 2., Cambridge 
University Press, 181 



Strauss, H.R. 1976, Pliys. Fluids, 19, 134 

Tang, H. Z. & Xu, K. 2000, J. Comp. Phys., 165, 69 

Toth, G. 1995, MNRAS, 274, 1002 

Troland, T. H. & Heilcs, C. 1986, ApJ, 339, 345 

Xu, K. 1999, J. Comp. Phys., 153, 334 

Zweibel, E. G. 2002, ApJ, 567, 962 



